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ABSTRACT 

Cold dark matter haloes are populated by high-density structures with sharply- 
peaked profiles known as caustics which have not yet been resolved by 3- 
dimensional numerical simulations. Here, we derive semi-analytic expressions 
for the density profiles near caustics in haloes which form by self-similar accre¬ 
tions of dark matter with infinitesimal velocity dispersion. A simple rescaling 
shows that similarly to the case of absolutely cold medium these profiles are 
universal: they are valid for all caustics and irrespective of physical param¬ 
eters of the halo. We derive the maximum density of the caustics and show 
that it depends on the velocity dispersion and the caustic location. We show 
that both the absolute and relative thickness of the caustic monotonically de¬ 
crease toward the center of the halo while the maximum density grows. This 
indicates that the radial component of the thermal velocities decreases in the 
inner streams, i.e. the collisionless medium cools down in the radial direction 
descending to the center of the halo. Finally, we demonstrate that there can be 
a significant contribution to the emission measure from dark matter particle 
annihilation in the caustics. 

Key words: dark matter, haloes of galaxies, caustics, velocity dispersion, 
dark matter detection 


1 INTRODUCTION 

Dark matter particles, if collisionless and cold, would focus under gravitational instability into caustics which are 
formally 2D manifolds of infinite density. In three-dimensional space, caustics are determined by tangent hyperplanes 
r = const to the phase surface in the six-dimentional phase space (r, v) and bound the regions of multi-stream flow 
where velocity has multiple values. Once formed, caustics of a given phase volume neither disappear nor overlap; a 
requirement of the Liouville theorem. However they can interact and merge with caustics of a different phase volume 
and consequently the nature of their singularity could change and they could undergo generic metamorphoses. In 
the case of potential flows, e.g. of light rays or cold dark matter on large scales, the singularities of the caustics 
and their metamorphosis have been classified up to three spatial dimensions (Arnol’d 1986, 1990). This classification 
remains intact in the presence of external or internal forces (e.g. in a self-gravitating system) for as long as the 
force is potential and smooth and the dark matter can be approximated as a collisionless fluid. The first and most 
common caustic has a density with an inverse square-root singularity (as occurs in the Zel’dovich approximation). 
This singularity has been rigorously proven to be robust in the case of a one-dimensional Vlassov-Poisson system 
(Roytvarf 1994) and is the only singularity that is of relevance to the present work. The other generic singularities can 
lie on one-dimensional manifolds (lines) or be isolated points. Generally, one-dimensional singularities are stronger 
than two-dimensional singularities and zero-dimensional are stronger than one-dimensional singularities. None of 
these singularities requires any particular symmetry for its formation (Arnol’d, Shandarin, Zel’dovich 1982, Arnol’d, 
Gusein-Zade, Varchenko 1985, Shandarin, Zel’dovich 1989). 

In cosmology, the study of formation and evolution of dark matter caustics has been historically and scientifi- 
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cally twofold: caustic related to the large-scale structure and caustic on the galaxy or smaller scales. On large scales, 
pioneering works by Zel’dovich and his collaborators (Zel’dovich 1970, Arnol’d, Shandarin and Zel’dovich 1982, Shan¬ 
darin and Zel’dovich 1989) showed that in a suitably-defined time and space coordinates, elements of a collisionless 
and self-gravitating fluid move on inertial trajectories: i.e. with their initial velocities. Thus, as in ray optics, when 
the path of the free-moving particles cross, density diverges, velocity becomes multivalued and caustics form. Caustics 
exist only as idealizations in models assuming that the medium is collisionless, continuous and cold i.e. the thermal 
velocity dispersion equals zero. It was shown that in the case of finite thermal velocity dispersion the density in the 
caustic regions becomes finite (Zel’dovich, Shandarin 1982). Although the real thermal velocity dispersion is never 
zero in many cases it is extremely small and the cold medium represent an excellent first approximation to reality. 
Thus, caustics of various types represent a very useful idealization for the study of complex density fields. It needles 
to say that both discreteness and collisionallity eliminate caustics. Caustics can bound regions with the morphology 
of filaments, sheets or clumps together forming a supercluster-void network which is remarkably similar to the mass 
distribution in the Universe on large scales (above 4 — 5 Mpc) as shown in numerical simulations ( e.g. Sahni and 
Shandarin 1996) and redshift galaxy catalogs {e.g. Bharadwaj et al. 2000). However, at small scales, within collapse 
structures, inertial approximation of Zel’dovich breaks down. 

The question arises whether caustics are relevant to the physics of dark matter at small scales, for example at 
the scale of a dark matter halo. Dark matter haloes can form for instance from the triaxial collapse of spherical 
perturbations or accretion of matter into the overdensities at the junctions of the filaments. In a cold dark matter 
Universe (with or without cosmological constant), they grow in a hierarchical manner by merging with other haloes 
and by accreting mass and are hosts to the formation of galaxies. In numerical simulations, they are bound overdense 
regions which are identified by various percolation algorithms. High resolution simulations, surprisingly, have found 
that the spherically-averaged equilibrium density profiles of cold dark matter haloes can be described by a power law 
(Dubinski & Carlberg 1991) and universal two-parameter function (Navarro, Frenk & White 1996, 1997). Although 
recent simulations have clearly established that many of the supposedly ’’relaxed” haloes still contain a large number 
of smaller subhaloes (Klypin et al. 1999, Moore et al. 1999) they do not show the presence of discrete flows and 
caustics (Moore 2001, Helmi, White & Springel 2003). We believe that they have not yet achieved enough mass 
resolution to observe the small-scale caustics in three dimensions and also probably suffer from spurious collisional 
effects (Melott et al. 1997, Splinter et al. 1998, Binney 2004) which wash out the caustics. 

It is worth remarking at this point that recent high resolution N-body simulations of the neutralino-dominated 
Universe have shown that a considerable number of the smallest haloes with masses as small as 10^ 6 Mq survive until 
the present epoch (Diemand, Moore & Stadel 2005). These simulations suggest that there must be about 10 15 such 
haloes in our galaxy. These haloes are expected to have very smooth caustics because there had been no smaller scale 
fluctuations in the initial spectrum. The results of our work are most directly applicable to this type of structures. 

Analytic evaluation of the halo density profile, and the prediction of the existence of caustics inside these struc¬ 
tures, started with the works of Gott (1975) and Gunn (1977) who used the spherically symmetric model. With the 
main objective of explaining the flattening of the rotation curves of the galaxies, they considered the formation of 
a dark matter halo from the secondary infall of matter onto an already formed galaxy (or in later works onto a 
spherical overdense region). In an Einstein-de Sitter Universe a spherical overdensity expands and then turns around 
to collapse. After collapse and at late times, the fluid motion becomes selfsimilar: its form remains unchanged when 
length are re-scaled in terms of the radius of the shell which is currently turning around and falling onto the galaxy. 
Physically selfsimilarity arises because gravity is scale-free and because mass shells outside the initial overdensity 
are also bound and turn around at successively later times. Self-similar solutions give power-law density profiles 
whose exact scaling properties depend on the central boundary conditions and on whether the fluid is collisionless or 
collisional (Fillmore and Goldreich 1984, Bertschinger 1985a, 1985b). The density profile obeys a power-law on the 
scale of the halo which provides an explanation of the flattening of the rotation curves of the galaxies. However, on 
smaller scales the density profile contains many spikes {i.e. caustics) of infinite density (with an artificial cutoff due 
to finite numerical resolution). The position and the time of formation of these caustics are among many properties 
which have been studied in the secondary infall model (Bertschinger 1985b). 1 

Such aforementioned studies not only have proven valuable for the prediction and description of large-scale 
structure, of dynamics and distribution of mass inside dark matter haloes and of galaxy formation but also recently 
for the detection of dark matter particles. Due to their significantly high density over their often already dense 
background, and their large number density, caustics are clearly of importance for dark matter search experiments. 

In the past few years, major experiments have gone underway for the direct and indirect detection of dark matter 


1 Although various elaborations have since been made on the secondary infall, in order to accommodate for the bi-scaling of 
the haloes density profile observed in simulations {e.g. see Henriksen 2004) here we concentrate on the original secondary infall 
model which yields a pure power law density profile. 
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particles. Direct detection experiments, such as DAMA and EDELWEISS, often use the annual modulation of the signal 
due to the orbital motion of the Earth around the Sun. Since the flux of dark matter in direct searches depends linearly 
on the local dark matter density, the search strategy and data analysis strongly depend on the spatial distribution of 
dark matter and its dynamics in the galactic halo. 

The indirect-detection experiments, such as ANTARES, HESS and GLAST, search for products of annihilation of 
dark matter candidates ( e.g. neutralinos) such as energetic neutrinos and 7 -rays. In the indirect searches, the flux 
of the annihilation products depends quadratically on the local dark matter density. Thus the degree of dumpiness, 
the density profile of a dark matter halo, the presence or absence of a central supermassive black hole and finally the 
presence of caustics, could influence the annihilation rate and boost the 7 -ray flux significantly. It has been shown 
that this boost is significant if there is a cusp at the centre of the halo (Stoehr 2003, Salati 2004). The accretion of 
dark matter into a central black hole if present in the halo, could also boost the gamma-ray flux by few orders of 
magnitude but again only if the dark matter profile develops a cusp at the centre (Gondolo and Silk 1999). Thus, a 
central core profile would not in general lead to a significant boost of the flux. However, although dark-matter-only 
simulations seem to show a cuspy profile in the centre of the haloes, some of the observations seem to contradict 
these predictions {e.g. see McGaugh et al. 2003). In addition, whether the cusp observed in the dark-matter-only 
simulations would survive in the presence of gas or would become less steep or else disappear due to the reaction 
from the baryonic gas, or is simply a numerical artifact (Binney 2004) is unclear. Caustics, on the other hand, would 
be inevitably present, as a direct consequence of Jeans-Vlassov-Poisson equation {e.g. see Alard & Colombi 2004 for 
a recent numerical simulations in 1 dimension). Therefore, it is worthwhile to study the density enhancement in the 
caustics and its possible implications for dark matter search experiments. Many properties of the dark matter halos 
related to detection of the dark matter signatures has been already discussed in early works. For example the velocity 
magnitudes of the peaks in the velocity space and the large-scale properties of galactic halos have been studied 
(Sikivie, Tkachev & Wang 1997) and the the geometry of caustics in the galaxy halos has been discussed (Sikivie 
1998, Sikivie 1999, Sikivie & Ipser 1998). A simple estimation of maximum density in caustics due to small thermal 
velocity dispersion has also been carried out (Bergstrom, Edsjo & Gunnarsson 2001). The formation and role of 
micropancakes in halos has been discussed (Hogan 2001). The more general question of the expected dimensionality 
of phase-space patterns in observations of galaxy structures has been studied and it has been suggested that the most 
prominent features will be stable singularities (Tremain 2001). 

The presence of a small velocity dispersion, e.g for neutralinos which are presently the most plausible dark matter 
candidates, smoothes the matter density at the caustic and gives it a finite maximum value. The principle problem 
addressed in this work is the precise derivation of the value of this quantity and its implications for dark matter 
search experiments. 

In view of the fact that almost all dark matter candidates have non-negligible velocity dispersion, we consider the 
secondary infall of dark matter with a very small but finite velocity dispersion. In this case the density profile at the 
very vicinity of the caustic would be affected and in addition caustics would have a physical cut-off to their density. 
Here, we evaluate analytic expressions for the density profile at the vicinity of the caustics and also determine the 
maximum density at the caustic positions. The analytic expression for the density profile of the caustics is given as 
a function of the initial velocity dispersion. The 7 - 1 'ay emission measure from the annihilation of the neutralinos in 
the caustics is then evaluated. Using our results, we evaluate the position, the thickness, the density and the 7 -ray 
emission measure for the first caustic of M31. 

This article is organised as follows. In Section 2, we review the basics of the secondary infall model. In section 3, 
we derive analytic expressions for the density profiles near the caustics in the presence of small velocity dispersion. 
I 11 Section 4, we use our density profile and evaluate a general analytic expression for the emission measure from a 
typical caustic. In Section 5 we evaluate the 7 -ray emission measure from dark matter annihilation in the first caustic 
(nearest to us) of M31. In Section 6 we conclude our main results. 


2 SELF-SIMILAR MODEL 

We consider a spherical overdensity of collisionless fluid in an Einstein-de Sitter Universe which eventually ceases 
expansion and turns around to collapse. The trajectory of a fluid element in radial motion obeys Newton’s law 

d 2 r Gm{r,t ) ^ 

dt? = 7 ’ W 

where the mass, m{r,t) inside a radius r, is not constant due to shell-crossing. At first the only way to tackle this 
problem seems to be via an N-body simulation. However, a major simplification arises once it is realised that the 
problem has a similarity solution [(Filmore & Goldreich 1984) and (Bcrtschinger 1985b) and we use the notations 
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of (Bertschinger 1985b) throughout]. The turn around radius rta(i) which is the only length scale in the problem is 
used to introduce the non-dimensional variables 


A = r( '\ ; e = lnf * ) 

?’ta (t) Vita/ 

; A/(A) = j m( © , 

4tt p H rl 

(2) 

where 



( t \ 8/9 

rta(t) = nta y -— J , 


(3) 


rita is the initial turnaround radius, tn a is the initial turnaround time, pH is the Einstein-de Sitter density {pH = 
1 /67T Gt 2 ) and fta is the turnaround time for a given particle ( i.e. when the particle is at its largest radius). In terms 
of the nondimensional variables (2) Newton’s equation (1) becomes 


d 2 A 7dX 

+ 9^ 


81 A = -9A^ 


(4) 


which has no explicit dependence on non-dimensional time £. The equation should be solved with the initial condition 
(at £ = 0 corresponding to t = tta) 


A = 1; 


d\ _ _8 

“ 9 ’ 


(5) 


and a prior knowledge of the mass M{ A). In the case of the Hubble flow, there is a simple solution to equation (4), 
before shell-crossing, which is given by M( A) = A 3 = Af ta e^ 2 ^ 3 . However, after shell-crossing occurs there are many 
particles having the same value of A. This can be taken into account simply by the summation 


M (A) = Afta^-l)* - ^ -24 * 73 

i 


(6) 


which adds (for i odd) the mass of the particles interior to A and subtracts (for even i) the mass exterior to A, 
accounting correctly for shell-crossing. Equations (4) and (6) can be solved numerically by iteration (see Appendix 
A for more detailed account). Here, we take a simpler approach. At small values of A (A < 1), mass becomes a 
power-law M{ A) « 11.2A 3 / 4 (Bertschinger 1985b). We take this fact into account and instead of solving (4) and (6) 
iteratively, use a simple approximation for M(A) 


M{ A) 


11.2A 3/4 
1 + A 3 / 4 ' 


(7) 


and then solve (4) numerically at the given gravitational potential generated by the mass distribution (7). 

As shown in Appendix A by Figure A1 the approximation (7) generates relatively small errors. A notable 
discrepancy between the approximation (7) and expression (6) appears only at the relatively large values of A ~ 1. 
However, for all the caustics under consideration the value of A is far less than one [the largest value of A for the first 
caustic is at A « 0.36 (see the table in Fig. Al)]. 

Solutions to equations (4) and (6), give a power-law density profile convolved with many sharp spikes (the 
caustics). As a particle expands to its turnaround radius, it collapses and re-expands again to its new maximum 
radius, which gives the time and position of the first caustics. It then re-collapses and re-expands to the position of 
the second caustic and so on (see Fig. A2). The calculation of the halo density profile itself is not the subject of this 
work. Here we are primarily concerned with the calculation of the density profile near the caustics and the maximum 
density at the caustics in the case of dark matter with finite velocity dispersion. 


3 MAXIMUM DENSITY IN COSMOLOGICAL CAUSTICS 
3.1 Cold medium 

First, we derive the equations of motion in terms of physical time t and radius r. From the definitions (2) one can 
easily obtain 


t = 


e ? , 

(8) 

r = 

V nta / 

Rta exp 

1 m 

(9) 

where 

Rta. = 

^ta(^ta) — ^*ita | 

<fta\ 8 / 9 

V £ita / 

(10) 
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is the turnaround radius reached by a particle at its turnaround time (which is equivalent to r( a in Bertschinger 
1985b). Solving equation (8) for £ one can obtain in terms of the physical coordinates and function A = A(£) the 
explicit solution 

/ t \ 8 / 9 f 

r(7? ta , t) = r ita y—J A i In 



Introducing dimensionless time r and dimensionless coordinates x and q 

t r R t z 

t = - , x = -, and q = - 

£ita ^ita ^*ita 

one can further simplify equation (11) 
x = r 8 ^ 9 A[ln(rq _9,/8 )]. 


( 12 ) 


(13) 


Equations (11) and (13) represent the mapping from Lagrangian space to Eulerian space parameterised by time. 
The density can be obtained from the conservation of mass 
i-i 


p(x) = p(q)- 


dx 

dq 


(14) 


where the ratio dx/dq must be taken at the time of formation of the caustic, Tk- The condition of the caustic formation 
dx/dq = 0 requires 

A' (&) = 0 (15) 

where A' = dX/dt) and = Tkq~ 9 ^ 8 denotes the position of a maximum of the function A(£). The derivative dx/dq 
at the Lagrangian distance A q from the caustic is 


dx 

dq 


d 2 x 

dq 2 


81 t 8/9 

Aq=—^ r \'/Aq, 
64 q 2 


(16) 


where A = A”(£fc) The relation between Aa; and A q can be easily found by expanding x(r,q) [equation (13)] into 
the Taylor series at the time Tk and using the condition of A = max [given by equation (15)] 


Aa; = 


A q + 


1 


d 2 x 
dq 2 


A 2 1 81 r 8/9 2 

= 2 64 ^ AfcA? - 


T k X 7 T k 

Thus, the inverse derivative (dx/dq)- 1 becomes 
i . „ / ,„\ -i/2 

(-Aa;)" 1 / 2 , 


dx 

dq 


1 8 q 


2 9 4/9 


2 


(17) 


(18) 


where the signs in the above equation reflect the signs of X'/ < 0 and Ax < 0 in the vicinity of the caustic. 
Substituting the derivative ( dx/dq ) _1 , obtained above, into equation (14) one derives the density in the vicinity of a 
caustic (additional factor 2 must be added due to two stream flow at Ax < 0) 


p(Ax) = A k (-Ax) 1/2 , 


(19) 


with 



( 20 ) 


where pH = l/6nGt 2 is the mean density of the universe, M ta = (3ir/8) 2 , Xk = X (£*,), and Xta = r ta (t)/rit a is 
the present dimensionless turnaround radius. [Equation (19) corresponds to equation (4.7) in Bertschinger (1985b).] 
Substituting r for x [equation (12)] one can obtain the density in terms of dimensional physical parameters. The 
parameters of the self-similar solution [£fc,-\(£fc) and A"(£*,)] that determine the density in the vicinity of every 
caustic must be obtained from numerical integration of equation (4) and equation (6) [or (4) and the approximate 
equation (7)]. 


3.2 Medium with thermal velocity dispersion 

In this section, we derive semi-analytic expressions for the density profile in the vicinity of the caustics for non-zero 
velocity dispersion following the method used in Zel’dovich and Shandarin (1982) and Kotok and Shandarin (1987). 
Although in rigorous mathematical terms caustics, defined as manifolds of infinite density, would not form in the 
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presence of a finite velocity dispersion, the density at the caustic “positions” would still be extremely high if the 
velocity dispersion is very small, as is the case for most dark matter candidates, and hence we still refer to these 
sharply-dense structures as caustics. 

The motion of a medium with small thermal velocity dispersion can be approximated as a simultaneous evolution 
of many streams with different initial velocities u at r = rta. The formation of the caustic in every stream occurs at 
different radii x v - We denote the distance from the caustic as 5x v = x — x v - We assume a linear relation between the 
relative position of the caustic and the initial velocity of the stream v 

Xv xo — akv, ( 21 ) 

where v is the dimensionless velocity related to physical velocity u as u = ( rit a /tita)v and ak is a negative constant 
to be determined numerically for every caustic. We will express the density as a function of the distance Air = x — xo 
from the caustic in the stream with zero initial velocity 

5x v = Ax — akV. (22) 


The major effect to consider in determining the maximum density in the ‘caustic’ in the medium with small thermal 
velocities is the shift of the position of the caustic in every stream with respect to the stream with v = 0. As a result 
the 1 /x /—Ax factor in equation (19) must be modified as in the following integral 

f{v)dv 


p(Ax) = 


(23) 


V —Ax + akv 

where f(v) is the velocity distribution function at the turn around radius q and the turn around time rta corresponding 
to a chosen caustic and f(v) is normalised to unity: J f(v)dv = 1 . The above integral is simply the sum of densities 
in all streams at a distance Ax from the true caustic in cold medium. 

In the simplest case of the one-dimensional top-hat (TH) velocity distribution 


/ th{v ) = | 2ct„ 


we obtain 


if |v| < <7„ 
otherwise 


(24) 


P( Ax) = 


Ak 


1 - 


- J- 1 - 


S\ 


Qt-kGv 


l - 


o 


Ax < -\a.k<Jv\, 
— \otk<Xv\ < Ax < | CXkCTv 

Ax > lafcerJ, 


(25) 


for the density profile near caustics, where Ak is given by equation (20). One recovers equation (19) from equation 
(25) in the limit of a v = 0. 

In order to find the constant ak we numerically solve the equation for A(£) in the presence of the initial nondi- 
mensional velocity perturbation 8\' Q in the initial conditions at £ = 0 so that 


A(0) = 1, 


A'(0) = — g+*Ai. 


The change of the maximum of A(£fc) = max = A k 
linear expression 

8Xk = AfcJAo, 


(26) 

Afc + SXk in this case is found to be well approximated by a 

(27) 


where the coefficient Afc depends on the caustic k and shall be determined numerically. Differentiating equation (13) 
one finds the relation 
dx 


dr 


-T - 1 / 9 


% m + a'(o 


(28) 


between the dimensionless velocity dx/dr and nondimensional functions. Therefore at the turnaround time, rta, the 
dimensionless velocity, v, given by equation ( 21 ), is 




5X' 0 . 


(29) 


Taking the variation of equation (13) at Tfc and using equation (27) one obtains the distance between the caustic in 
the stream with velocity v and the caustic in the stream with v = 0 

x v - xo = T% /9 A k 5X' 0 . (30) 

Finally, recalling that x v — *o = akV [equation (21)] one obtains the expression for afc in terms of Afc 
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otk « T ta exp Afc, (31) 

noting that Afc < 0 and ak < 0. Combining equations (20),(25) and (31) gives the full expression for the density in 
the vicinity of caustics for the top-hat velocity distribution function. We derive the density profile in the vicinity of 
a caustic for the exponential and Gaussian velocity distribution functions in Appendix B. 

A much simpler and approximate expression for the density can be written as 



' (—Ax)“ 1/2 

for Ax < Ax c k, 


p{ Ax) = Ak < 

(-Ax c fc)- 1/2 

for Ax c k < Ax < 0, 

(32) 

where 

• 0 

for Ax > 0. 


== Otj^CTv 

f 8 c ' 

= (T v r ta exp 

V9 / 

) Afc, 

(33) 


note that A x c k is negative. The above expression approximates the maximum density for the exponential and Gaussian 
velocity distribution functions quite well (see Fig. 1). Physically this approximation means that the highest density 
in caustics is reduced by thermal velocities but since the amount of mass having the highest densities is small 
oc (—A Xck) 1 / 2 oc it does not significantly affect the rest of the density distribution. Although a similar idea was 
used in the estimate by Bergstrom, Edsjo, Gunnarsson (2001) they have not derived the thickness of the caustic (eq. 
33 and 37) and have not shown the accuracy of their approximation. 

In order to evaluate the density profile in the vicinity of a caustic we numerically solve equation (4) with the 
mass distribution approximated by equation (7) and the initial conditions given by equation (26). We summarise the 
caustic parameters (including Afc) in Table 1. 


k 


Afc = A (£fc) 

A l = A "(&) 

Ak 

1 

0.985 

0.368 

-5.68 

-0.070 

2 

1.46 

0.237 

-11.2 

-0.025 

3 

1.76 

0.179 

-16.7 

-0.014 

4 

1.98 

0.146 

-22.3 

-0.0085 

5 

2.16 

0.124 

-28.0 

-0.0059 

6 

2.31 

0.108 

-33.9 

-0.0043 

7 

2.43 

0.0960 

-39.8 

-0.0034 

8 

2.55 

0.0866 

-45.7 

-0.0027 

9 

2.64 

0.0790 

-51.7 

-0.0022 

10 

2.73 

0.0728 

-58.8 

-0.0018 


Table 1. The nondimensional parameters of the first ten caustics obtained from the numerical fits to the curves in the phase 
space (similar to those shown in Fig. A3) are shown. 

Rearranging (32), we arrive at the following full expression for the maximum density (which occurs at Ax = Ax c fc) 


Pmax — 


V ^ tt 2 e - 11 ^/ 18 1 

8 y/Apr 



(34) 


where t is the present time. Factor (t/tita) -1 ^ 18 ^ 1 ^ 2 can be expressed in terms of radial component of the physical 
thermal velocity dispersion, a p h, at the turnaround radius Rta at the turnaround time ft a - First, we note that 
Ov = <fph(nta/tita) _1 , and then we use equation (3) to obtain nta/hta = ita) 1 ^ 9 - Combining both factors 

one obtains 



We wish to stress that both t and rta are present time and present turnaround radius while <r p h is the physical 
velocity dispersion at the turnaround radius at the turnaround time. The thermal velocity dispersion at turnaround 
radius at the turnaround time can be estimated from the conservation of the phase space volume. The density at the 
turnaround radius is D( 1) = (37r/8) 2 « 1.39 times greater than prr(fta) therefore <r p h approximately D(l) 1 ^ 3 « 1.24 
times greater then the thermal velocity dispersion in the homogeneous universe at that time. Thus, the maximum 
density only depends on the caustic, the velocity dispersion and the background Einstein-de Sitter density, pH- The 
ratio of the maximum density to the background (Einstein-de Sitter density which should not be confused with the 
local halo density) is thus almost independent of any physical parameters. Using caustic parameters given in the 
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e 

Figure 1. The rescaled density r/ = pA^ 1 [c\p(7 v ) 1 versus the rescaled distance e = AxaZ 1 ^ 1 from the caustic. Dashed line 
shows the case of cold matter (a v = 0), circles show the top-hat velocity distribution. Dashed-dotted and solid lines show the 
cases of the exponential and the Gaussian initial velocity distributions respectively. The dotted curve is a simple approximation 
given by equation (32), which is used in the evaluation of the emission measure in Section 4 and used also in Figure 3. The 
density profile is evidently universal: it is independent of halo parameters and is valid for all of the caustics. 


previous table, we can evaluate the maximum density at caustic positions. In Table 2, we summarise the value of 
p/ Ph / y/ov and also of the local halo density for the first ten caustics. 


k 

i 

2 

3 

4 

5 

6 

7 

8 

9 

10 

pmax/ pH / y/d v 

li. 

24. 

39. 

56. 

74. 

95. 

117. 

139. 

165. 

190. 

Phalo/ PH 

12. 

40. 

83. 

139. 

210. 

297. 

397. 

506. 

641. 

777. 


Table 2. This tabic gives the maximum caustic density evaluated using (32) at the first ten caustic and also the halo density 
evaluated using the approximate expression (7) at the position of the caustics. Both of these densities are given as a ratio to 
the background Einstein de-Sitter density, pp . The non-dimensional velocity dispersion, a v , is given by expression (35). 

Thus, we can evaluate for a given velocity dispersion, the positions when the maximum caustic density becomes 
equivalent to the background density. Clearly for small values of a v this would occur only at small radii for inner 
caustics and vice versa . The values of velocity dispersion for cold dark matter is very small and would be expected 
to be much smaller than unity. Thus, it is clear from Table 2 that the enhancement factor can be extremely high, for 
low values of the velocity dispersion a v . 

To summarise this section, we also write the approximate maximum density of the caustics (see the “Approxi¬ 
mate” profile of Fig. 1) and their thickness obtained by our method. We have for the maximum density of caustics 


Pmax — 


/ ^V3^g-17W18^ 


rta(f) 


pH it). 


\y2V2 3 1 / 3 y'AfcA l K) V tcr ph (t) 

The radius of the caustic shell and its thickness in physical units are 


rk = Afcrta(t), Ar fc = —e S5fe/9 A k ta ph (t) 


(36) 


(37) 


It is remarkable that the thickness of the caustics is universal and depends only on the present day DM velocity 
dispersion and not on the mass of the DM halo. 
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Figure 2. Main plot: The top line (filled circles) shows the separation of caustics, A k — Afc_|_i, as a function of the distance, 
A*,, from the center. The two bottom lines show the distance between caustics in the streams with negative and positive 
perturbations of the initial velocity (<5 Aq in eq. 26). Note the logarithmic axes. The two lower lines correspond to <5 Aq = ±0.001 
(crosses) and <5 Aq = ±0.01 (open squares). All three lines can be well fitted by a power law oc A?' 1 showed by the dashed 
line. Inset: The ratio of the thicknesses of the caustics to their separations , 2<5Afc/(A; c — Aj,_|_i), is shown as a function of their 
radii. The vertical axis is linear. The ratios converge to approximately 0.0007 and 0.0066 for the two different values of velocity 
dispersions as marked on the plots, which shows its linear dependence on the initial velocity perturbation. 


Both the absolute and relative thickness (Ar*, and A rk/rk) of the caustic monotonically decrease toward the 
center of the halo. This indicates that the radial component of the thermal velocity decreases toward the center. We 
do not consider the evolution of the angular components of the thermal velocity in this paper, however we would like 
to speculate that they grow toward the center making the velocity distribution function anisotropic (oblate ellipsoid). 
The decrease of the radial temperature can also be viewed as a consequence of the Liouville theorem that forbids the 
overlapping of the streams in the phase space. 

Figure 2 shows the logarithm of separation of the caustics, A*, — Xk+i, and thicknesses, 25Xk, [see eq. (27)] and the 
ratio of these two quantities (plotted in the inset) as a function of the logarithm of the radius, \k . Both the separation 
and thickness are scaled as a power law oc A 2 ' 1 : the formal fit to the separation of caustics is Xk — A*,+i = 1.24A 21 and 
for the two examples in the main plot (corresponding to the initial velocity perturbations SX' 0 = ±0.01 and ±0.001) 
this fit is re-adjusted by factors of 0.066 and 0.007 respectively (see the inset). This is a remarkable result which shows 
that in the course of gravitational evolution the streams remain well-isolated from each other in spite of the fact that 
their separations diminish. 


4 EMISSION MEASURE (£M) FROM DARK MATTER PARTICLE ANNIHILATION IN THE 
CAUSTICS 


The annihilation flux (in photons/cm 2 /s) can be written as 


x± / dn x (sm) 


Afl 


( 38 ) 
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caustic surface 



Figure 3. An approximation is used in which the caustic is the outer-shell (larger circle) outside which the density is zero, 
inside a layer of thickness Ax c k the density is constant taken to be p c = A^Ax ck ' . In the interior of this shell the density 
falls as l/V —Ax where Ax is the radial coordinate measured from the tangent point D (Ax < 0). 


where the emission measure, 

EM. = / p 2 (s)ds , (39) 

J line of sight 

is found by integrating the square of the density along the line of sight and over the solid angle Afl, m 7 is the mass 
of the candidate particle ( e.g. neutralino), N 7 is the number of photons produced per annihilation. To compute the 
first part of integral IV 7 (av) / (Aimi^) a supersymmetric model needs be selected. We shall not discuss this aspect 
here which is already extensively discussed in the relevant literatures (e.g. see Jungman, Kamionkowski and Griest 
1996 and references therein). In this article, we obtain analytic expression for the emission measure (39), using our 
approximate expression for the density (32). 

Here we are interested in calculating the (£M) from the regions close to the caustic surfaces formed in cold 
matter. Clearly, the (£M) is considerably higher when the line of sight is close to the tangent to the caustic surface. 
We estimate the (£M) in a small vicinity of this tangent. Although, in principle in order to obtain the emission 
measure (39) one can integrate expressions (B6) or (B7) or (B8) numerically, here we use our simple approximation 
(32) for the density profile in the vicinity of the caustic and make analytic estimates for the boost factor (39). 

We assume that the density in the vicinity of the caustic can be approximated by equation (32). Figure 3 
illustrates the geometry of the system. The figure shows the plane passing through the observer, O , the centre of 
the galaxy, C, and point D where the line of sight OP is tangential to the caustic sphere. The external and internal 
caustic spheres have the radius R ex and I?j n = 7? ex + A x c k respectively (note Ax c k is negative, see equation (?? and 
also Fig. 3). The density in the shell between two spheres is constant p = Ak(— A* c fc) -1 ^ 2 , while inside it falls as 
p = Ak(—Ax) -1 / 2 where Ax is the radial coordinate measured from point D on the external sphere (Ax < 0). 

We evaluate the (£M) as a function of the angle 6 measured from the line OP upwards. First, we calculate the 
contribution to the emission measure along the line of sight which runs inside the angle POQ. In this case, no integral 
needs to be evaluated since the density is constant and the integral of density along the line of sight is the density 
times the length of the chord between two points where the line of sight crosses the external circle (e.g. tT). 

The contribution to the emission measure, (39), from the lines crossing inside POQ is 


£M(0) = 


p 2 ds = p 2 c L ch (6), 


(40) 


where L c h(9) is the length of the chord. The equation of the circle in the vicinity of the tangential point D is to the 
lowest order 

,2 ; _ 

Ax = — or y = ±\Z~2R^Ax (41) 

^JXex 

where y is the coordinate along DP measured from point D. For example, the length of the chord tT is 

Ltr = \J(yr + yt ) 2 + (Axt — A xt ) 2 ~ yr + yt, (42) 
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because Ax is of the higher order in y and thus the second term can be neglected (yr > 0 and yt < 0). Neglecting 
the difference between A xt and A xt allows to approximate them as A xt ~ A xt ~ —Lod9 = —Yd9 (where 
Y d = sjY 2 - XL is the distance from the observer to the caustic and Yc is the distance from the observer O to the 
centre C) assuming 9 is small. Combining this approximation with the second equation (41) the length of the chord 
can be written in a simple form 

L ch = 2V2 RexYoVd, (43) 


and thus, 

£M{9) = 2plV2R ex Y D V9 for 9 <9 C = LPOQk. ~^ Xck . (44) 

Yd 

Next, we suppose that the line of sight is along OS in which case it crosses the region where the density falls as 
(—Ax) -1 / 2 . The integral (39) is now the sum of three parts: 

SMhH = (£A i)hb + ( £M)bh + ( £M)bB = (£M)i ,b + pc(Lhb + Lbh). (45) 


The length Lhb + Lbh = Lhh — LbB can be easily evaluated in a similar manner as before. The length of the chord 
hH is given by equation (43) while for the chord bB one needs to substitute R\ n = R ex + A x c k for R ex , \JYq — Rf n 
for Yd, and the angle <j> = LQOS for 9. For small angles 9 and small Ax c t, the change in the radius and distance to 
the caustic result in the higher order corrections and can be neglected, yielding 

LbB = 2 \J 2(7iex + Ax c k)\/Y£ — R'fi 1 \f4> « 2V2R ex Y D ^9 — 9 C . (46) 

where 9 C = — Axcu/Yd = Afi./ (p^Yo) and all distances are given in units of equation (13). Thus, the (£A4) from the 
parts of the line of sight with constant density becomes 

(£M) hb + ( £M)bh = 2p 2 c V2R ex Y D (VO - yJJG ofij for d> 9 c = -^£. (47) 

The complicated part of contribution to the integral (39) comes from the line of sight bB. Here, we actually need 
to calculate the integral (39) where the density is no longer constant but falls as (—Ax) -1 / 2 . First of all we need to 
express the distance |Ax| = L e p in terms of the line of sight distance a = Lhf- This can be easily done by solving 
the intersecting chords relation LhF x Lfh = L e F x Lfe for Ax (Lfe = 2 R ex — |Ax| and Lfh = L — z (where 
L = L^h) which yield 


ZB 

2 f 


( £M)bB = pc p(z) dz = p c 


dz 


Re X VRle X LbBZ T Z^ 


where Zb and zb correspond to points b and B respectively. The integral (48) can be written in close form 


I = I—-± ^ arctanh 

jL" 


2 R ex L(L — 2z)VRix — Lz + z 2 


-In 


z(L - z){AR 2 ex + L 2 ) - 2R 2 X L 2 ) 

The limits in the integral are to linear order 

= V2R ex Y D ^ \fi§ — \J9 — 9^j and zb = \Z2R ex Yu ^Vd + \J9 — 9^j . 


| — hi ^2\/R 2 X — Lz + z 2 — L + 2z^ 


Zb 


(48) 


(49) 


(50) 


Substituting the limits into equation (49) is straightforward but results can be lengthy and complex. In addition, 
despite the exact form of the integral (49) we use the accurate limits only to the lowest order. Thus, we simplify the 
equation for (£M) by taking the lowest order terms in the series expansion for small 9 and 9 C . 


(£M)(9) = p c 


2R ex i f \[9 + V 9 — 9 C 


: In 


Yd fi9 \V9 - V9 ~T c 
Thus, collecting various expressions for the emission measure together, we arrive at 


£M= pi 


0 for 9 < 0, 

2\/2 Re X YoVd for 0 < 9 < 9 C 
e-6 c 


2V2Re X Yo (Ve - V9~efi + yj2R ex /Y D 9~ 1/2 In for 9 > 9 C . 


(51) 


(52) 


where once again we mention that R ex is the radius of the caustic sphere, Yd is the distance to the caustic surface, 
9 C = —Axck/Yo = A^Kp^Yd) and all distances are given in units of equation (13). 
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Figure 4. Contribution to the emission measure (39) along different lines of sights as marked on Fig. 3. All values are in 
non-dimensional coordinates. 


5 EMISSION FROM A CHARACTERISTIC M31-TYPE DARK MATTER HALO 

Detailed application of our results to haloes especially minihaloes (Diemand, Moore & Stadel 2005) shall be presented 
in a forthcoming work. Here, as an exercise, we apply the results of the previous sections to a very simple model of 
M31. Emission from M31 dark matter halo has already been studied, largely in the context of dark matter search 
experiment CELEST ( e.g. see Nuss et al. 2002 and Falvard et al. 2004). M31 is situated at a distance of 780 Kpc from 
the centre of MilkyWay (MW) with a present turnaround radius which is taken to be at about 800 Kpc (Sandage 
1986, Karachentsev 2002). In fact M31 and MW can be considered to be embedded in a common halo, however here 
in our very simple model we consider M31 to have its own halo. The first caustic, which is the closest to us cuts the 
radius joining the centre of MW to the centre of M31 at a distance of about 500 Kpc from us. The tangent point 
to this caustic lies at about 700 Kpc from the centre of MW. This caustic has a thickness of about 0.115 Kpc and 
subtends an angle ( 9 C ) of about 0.006° at the centre of Milky Way. The maximum emission measure from this caustic, 
using the expressions (52) for a velocity dispersion of a v = 0.001 is given in Table 3. We assume a Held of view of 1°, 
which is at the lowest-end for most detectors. 


k 

r (kpc) 

Sr (kpc) 

0 C ° 

Pm ax 

C c £Mde 

lEM dO 

Total (EM) Gev 2 cm 5 c 4 

~T 

300 

0.115 

0.006 

350pfr 

260plfrta 

5 x 10 5 p 2 ff rta 

4 x 10 24 


Table 3. This tabic gives approximate values for various caustic parameters and also the emission measure for the first outer 
caustic of M31. which is the nearest to Milky Way. We take the non-dimensional velocity dispersion to be a v ~ 0.001 which 
would give only a very modest estimate of the maximum density. For neutralino, this parameter is smaller by a few orders of 
magnitude leading to a significantly higher density maximum but at the same time smaller angle, 9 C . Expression (35) can be 
used to transform between the physical velocity dispersion, o p h, and the non-dimensional velocity dispersion, a v . 

The emission measure from the first caustic of M31 is at least of the same order 2 as that for example from the 

2 Different values for the emission measure form the galactic centre has been evaluated ranging from 10 21 to 10 31 depending on 
various physical assumptions such as the presence or absence of a central core, a central cusp, or a central supermassive black 
hole (see for example Stoehr et al. 2004; Evans, Ferrer and Sarkar 2004). 
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centre of Milky Way. Here, we have used a moderate value of a v = 0.001, a lower value would sharply increase the 
maximum density as is evident from expressions (32) and (38) although it would also reduce the thickness of the 
caustic. The simple exercise in this section also demonstrates that emission measure from caustics can serve as a 
mean to put bound on the mass of the dark matter particle candidates. 


6 DISCUSSION 

In this work, we have evaluated the density profile near caustics which arise in the selfsimilar scenario of the formation 
of dark matter haloes. We have obtained a universal analytic expression for this density profile and its maximum 
value in the presence of a small velocity dispersion. We have shown that the maximum density at the caustics depends 
primarily on the velocity dispersion and is significantly higher than the halo density at the position of the caustic 
for outer caustics. The radius at which the caustic density approaches the background halo density depends on the 
value of the velocity dispersion and is expected to arise only for the innermost caustics. We have then evaluated 
the emission measure from the caustics and applied it to the concrete example of M31 type halo outermost caustic, 
which is closest to us. This example demonstrates that caustics can be promising sources for dark matter search 
experiments. Application of our results to other haloes, in particular to small haloes which have not had significant 
mergers, to specific search experiments such as HESS and also to direct detection experiments remains to be done. 

Although we have considered an Einstein-de Sitter Universe, we expect our results to give reasonable approxi¬ 
mation for the ACDM Universe as well. The role of dark energy becomes significant at rather small redshifts (~ 0.3) 
which we expect to be well after the formation of the typical dark matter haloes we consider here. Furthermore, 
once a particle turns around and collapses, it separates from the background expansion and its subsequent motion 
should not be effected by the A term. However, one ought to use real density of dark matter halo and not the critical 
values. The second caveat in our consideration is the assumption of spherical symmetry of haloes and their cold 
accretion which does not hold in general cases. However, even in real 3D collapses, parts of the caustic surface can be 
well-approximated by spheres. Finally, the third and probably most serious problem is associated with the smooth 
precollapse conditions on galactic scale that contradicts the hierarchical clustering scenario. The cold dark matter 
models predict a relatively high level of small scale perturbations that result in the formation of small gravitationally 
bound haloes that are assembled into more massive haloes at later times. Thus, the dark matter accrets onto halos of 
galactic size in the form of smaller haloes that may significantly affect the density in the vicinity of caustics. However, 
there is a possibility that in cold dark matter models the smallest haloes can survive the tidal destruction in more 
massive halos (Diemand, Moore & Stadel 2005). In this case our results can be applied directly. This issue will be 
studied in detail in the following work. 

We show that both the separation of neighboring caustics and their effective thicknesses scale as a power law of the 
radius oc A^' 1 . This scaling demonstrates that the streams in the phase space corresponding to different macroscopic 
velocities remain well isolated despite of the fact that the velocities of inner streams vanish at the center of the halo. 
The radial component of the microscopic thermal velocities also vanishes as the stream descends to the center. We 
call this effect gravitational cooling of the radial temperature. Our method can be generalized for the case of nonradial 
components of the thermal velocity and we address this question in the following work. 

We believe that the current results represent a step toward building a more comprehensive theoretical model of 
the gravitational collapse. A more elaborate study of caustic distributions and density profiles, in broader settings 
without the assumptions of spherical symmetry, inertial trajectories or smooth initial conditions remains a challenging 
task. 
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APPENDIX A: NUMERICAL SIMULATIONS 

There is a computational difficulty in solving equation (4) which is singular at A = 0 where the velocity becomes 
infinite and changes discontinuously from negative to positive. In order to avoid this problem a small amount of 
angular momentum, of the form J 2 /A 3 is added the RHS of equation (4) where hereafter the value of J = 1CP 9 is 
adopted (Bertschinger 1985b). 

Further computational issues arise in solving (4) and (6). The outline of the numerical procedure is as follows: 
a guess is first made for the mass distribution e.g., M(X) = Mt a A 3 ' 4 for A < 1. Then equation (4) is integrated to 
obtain A(£). This is used to obtain a new approximation to M( A) and so on until a self-consistent solution is found 
(see Bertschinger 1985b for full details). In this article we have chosen the fitting formula (7) instead of (6), which as 
shown in Fig. Al works well for small values of A. 

The solutions to equations (4) and (7) is plotted in Fig. A2 which can be viewed as the trajectory and phase 
diagram of one particle during the course of evolution of the halo, or as a snapshot of the positions of many particles 
in the halo. The change in real and phase space for finite velocity dispersion is clear: velocity dispersion leads to 
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k 

5 


d\/d £ 

1 

0.985 (0.988) 

0.368 (0.364) 

-5.86 (-6.35) 

2 

1.46(1.46) 

0.237 (0.236) 

-11.2 (-11.6) 

3 

1.76(1.76) 

0.179 (0.179) 

-16.7 (-16.9) 

4 

1.98 (1.98) 

0.146 (0.145) 

-22.3 (-22.4) 

5 

2.16(2.15) 

0.124(0.123) 

-28.0 (-28.1) 

6 

2.31 (2.30) 

0.108 (0.108) 

-33.9 (-33.8) 

7 

2.43 (2.43) 

0.096 (0.096) 

-39.8 (-39.6) 

8 

2.55 (2.54) 

0.087 (0.086) 

-45.7 (-45.6) 

9 

2.64 (2.64) 

0.079 (0.079) 

-51.7 (-51.5) 

10 

2.73 (2.73) 

0.073 (0.073) 

-57.8 (-57.6) 



Figure Al. The fit (solid line) M(A) = 11. 2A 3 '/ 4 ' /(l. + A 0 ' 75 ) to the data points (open circles) taken from Table 4 of 
Bertschinger 1985b. The above table shows our values and Bertschinger 1985b values (given inside brackets) for A, £ and 
d 2 X/d£ 2 for the first ten caustics. 


the broadening of caustics. The density at caustics no longer diverges but has a maximum cut-off determined by the 
velocity dispersion, which is the main result of this work. 

The small velocity dispersion effects are primarily related to defocusing of trajectories of particles with different 
thermal velocities in the vicinity of caustics. Imagine the evolution of large number of streams each corresponding 
to a particular value of thermal velocity then these streams will produce a caustic at a slightly different radii. The 
resulting density field becomes the sum of densities in every stream. We assume that each stream evolves in the same 
gravitational field generated by the mass distribution of the cold medium [equation (7)]. 

We consider the process in the nondimensional coordinates £, A, X'dX/d £. As we have mentioned in Subsection 
3.2, the effects of thermal velocity dispersion can be considered by adding small velocities <5Ag = dA^O) to the initial 
velocity. Thus, the initial conditions become 

Ao = A(0) = 1, Aq = A (0) = —— + SXq. (Al) 

The result of the integration is shown in Fig. A3 for the first caustic. For small values of 5X' 0 , the major effect on the 
caustic is the change of the maximum value Xk that can be well-approximated by a linear function 

SXk — A k 6X 0 (A2) 

as is evident from Fig. A4. We find A*, by fitting the numerical results for every caustic k = 1,..., 10. 


APPENDIX B: THE DENSITY PROFILES FOR THE TOP-HAT, EXPONENTIAL AND 
GAUSSIAN VELOCITY DISTRIBUTIONS 


We evaluated the density profile in the vicinity of a caustic for the top-hat initial thermal velocity distribution 

/th(u)=po(^ \v\<a v ( B1 ) 

10 otherwise 

in Sec. 3.2 [equation (25)]. Using similar approach one also can derive the density profiles for the exponential and 
Gaussian velocity distributions 


Ie(v) 


1 

Pot ,— ex P 
Z(T V 



(B2) 
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Figure A2. The left plot is the non-dimensional trajectory A(r) for a collisionless fluid, with minute velocity dispersion, given 
by similarity solution to equations (4) and (7). The particles reach their maximum radius at A = 1. The non-dimensional 
coordinate A and the non-dimensional time £ at the maxima depend almost linearly on the initial velocity of the particle. The 
right plot is a small portion of the phase diagram given by the similarity solution. Each particle travels along the entire curve 
and at a given time there is a particle at each point on the curve. As we go to smaller and smaller radii, the strips become 
narrower and more closely packed and finally resemble a smooth distribution. To demonstrate the thickening of the caustics we 
have taken an unrealistically large range of velocity dispersion, which is causing the shift, seen in the left panel of this figure, 
in the formation times of the caustics. 


fa{v) = po 


\J 2 7T (7 v 


exp - 


2 <rjj 


(B3) 


Thus, deriving the density profile in the vicinity of a caustic in both cases consists in straightforward evaluations 
of a few integrals. However, one can simplify the calculations by introducing the scaled distance from the caustic e 
and the scaled density 77 

Ax 


5 V ~~ - 

\ a k<Jv\ Ak\akav \~ 1 '' 2 

In terms of these variables the density of cold dark matter with zero velocity dispersion has a simple form 
f (—e ) -1 / 2 for e < 0 


r)a v = 0 


(B4) 


(B5) 


0 


for e > 0 


that can be obtained as a limit at cr„ —> 0 of any of three expressions bellow. 

For the top-hat velocity distribution function, with non-vanishing u v , we have 

for e < — 1 , 
for - 1 < e < 1 , 

0 for £ > 1 . 

For the exponential velocity distribution function one obtains 

{e _£ [l — erf(yr^)j _ j e £ erf(i-s/T7)} f or e < 0 



(B 6 ) 
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vfzio- e 


(B7) 


for e > 0 


and finally for the Gaussian velocity distribution function the density is as follows 
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Figure A3. The flipped plot of the phase space around the first caustic for three values of the initial velocity perturbation, 
<5 Aq. Phase plots are used to evaluate the density at the caustics in the presence of a small velocity dispersion. 





Figure A4. The left plot shows the variation of the non-dimensional position of the first caustic with the initial velocity 
perturbation, ^A'fO). The middle and right panels show the variations of the second derivative and the time of formation of 
the first caustic with velocity dispersion, respectively. The quadratic fits (solid line) are also shown. The dashed straight lines 
demonstrate that for very small velocity dispersion, a linear function would be equally appropriate. Similar fits were obtained 
for the first ten caustics. 
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vg = 


e / 4 27rBesselI j, — ^2BesselI< 
^X|e- 2 /4 B esselK(i,4) 


for e < 0 
for e > 0 . 


The density profiles r) = 77 (e) are shown in Fig. f. 

We also show a simple approximation used in further calculations 

(- e )- 1/2 


for e > 0 . 
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